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Abstract 

We propose a framework for studying predictability of extreme events in complex systems. 
Major conceptual elements — hierarchical structure, spatial dynamics, and external driving - 
are combined in a classical branching diffusion with immigration. New elements — observation 
space and observed events — are introduced in order to formulate a prediction problem patterned 
after the geophysical and environmental applications. The problem consists of estimating the 
likelihood of occurrence of an extreme event given the observations of smaller events while the 
complete internal dynamics of the system is unknown. We look for premonitory patterns that 
emerge as an extreme event approaches; those patterns are deviations from the long-term sys- 
tem's averages. We have found a single control parameter that governs multiple spatio-temporal 
premonitory patterns. For that purpose, we derive i) complete analytic description of time- and 
space-dependent size distribution of particles generated by a single immigrant; ii) the steady-state 
moments that correspond to multiple immigrants; and iii) size- and space-based asymptotic for 
the particle size distribution. Our results suggest a mechanism for universal premonitory patterns 
and provide a natural framework for their theoretical and empirical study. 
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I. INTRODUCTION 



Extreme events are a most important yet least understood feature of natural and socioe- 
conomic complex systems. In different contexts these events are also called critical tran- 
sitions, disasters, catastrophes, or crises. Among examples are destructive earthquakes, 
El-Ninos, heat waves, electric power blackouts, economic recessions, stock-market crashes, 
pandemics, armed conflicts, and terrorism surges. Extreme events are rare, but consequen- 
tial: they inflict a lion's share of the damage to population, economy, and environment. 
The present study is focused on predicting individual extreme events. That problem is piv- 
otal both for fundamental understanding of complex systems and for disaster preparedness 
(see e.g., 

Our approach to prediction is complementary to more traditional and well-developed 
ones, which include classical Kolmogoroff- Wiener extrapolation of time series 4], 5|, linear 
(Kalman-Bucy) |6j and non-linear (Kushner-Zakaf , [lo| filtering, sequential Monte- 
Carlo methods [9], or the extreme- value theory The need for a novel approach is 
dictated by a non-standard formulation of the prediction problem, where one is partic- 
ularly interested in the future occurrence times of rare events rather than the complete 
unobserved state of the system in continuous time. We notice, accordingly, that often the 
easily observed extreme events can not be defined as the instants of threshold exceedance 
by the observed physical or economical fields, like air temperature or asset price. A paradig- 
matic example is an earthquake initiation time, which is determined by complex interplay 
of stress and strength fields in the heterogeneous Earth lithosphere. The physical theory 
for spatio-temporal evolution of these fields is still in its infancy, their values can hardly be 
measured with the existing instruments, or predicted using the available statistical meth- 
ods. At the same time, earthquakes are readily defined, measured, and studied. 

Prediction here is based on analysis of observable permanent background activity of the 
complex system. We look for premonitory patterns, i.e., particular deviations from long- 
term averages that emerge more frequently as an extreme event approaches. These patterns 
might be either perpetrators contributing to triggering an extreme event, or witnesses merely 
signaling that the system became unstable, ripe for a disaster. An example of a witness is 
proverbial "straws in the wind" preceding a hurricane. 

The following types of premonitory patterns have been established by exploratory data 
analysis and numerical modeling: (i) increase of background activity; (ii) deviations from 
self-similarity: change of the size distribution of events in favor of relatively strong yet 
sub-extreme events; (iii) increase of event's clustering; and (iv) emergence of long-range 

ogy and 
Im- 



correlations. Solid empirical evidence for existence of these patterns in seismo 
other forms of multiple fracturing has been accumulated since the 1970s 0, Q 
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portantly, these patterns are universal, common for complex systems of distinctly different 
origin. Similar premonitory patterns have been observed in socio-economic systems 39j,|40j], 



dynamic clustering in elastic billiards 42], hydrodynamics, and hierarchical models of ex- 
treme event development 43 48|. We propose here a general mechanism that reproduces 



these universal premonitory patterns. 

We focus in particular on premonitory deviations from self-similarity. Self-similarity 
is one of the most prominent features of complex systems. A canonical example is a 
power-law (self-similar) distribution of system's observables, whose remarkable feature is 
inevitability of extremely large events that dwarf numerous smaller events. Power-law 
distribution is well known under different names in such diver se p henomena as inertial- 



range self-similarity in turbulence (Kolmogprov-Obukhov laws) |49l453|. energy released in 
an earthquake (Gutenberg- Richter law) [54j-l56| , word usag e frequency in a language (Zipf 
law) |57j, allocation of wealth in a society (Pareto law) |58j,|59|, war casualties (Richardson 
law) 
slide 
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number of papers published by a given scientist (Lotka law) 16111 . mass of a land- 
number of species per genus [671 ]. and many other 



63J, stock price returns 
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An important paradigm of self-organized criticality |73j, |74j that is demonstrated 
by sand-pile 75[, forest-fire [76J, and slider-block 77H79j models and their numerous ramifi- 
cations has been introduced in order to understand dynamic processes whose only attractor 
corresponds to self-similarity (criticality) of the size distribution of appropriately defined 
events. 

Exact self-similarity, as well as many other universal properties, however is only an 
approximation to (or a mean-field property of) the observed and modeled systems; at each 
particular time moment the distribution of event sizes deviates from a pure power-law 
form. We show in this paper how to use such deviations for understanding the dynamics 
of a complex system in general and occurrence of extreme events in particular. 

The rest of the paper is organized as follows. We informally outline our model and 
the corresponding prediction problem in Sect. [Til A formal model description is given 
in Sect. II III Section [IV] summarizes the study's results most relevant to the prediction 
problem. Section [V] derives the spatio-temporal model distribution as a function of the 
control parameter. Section [VTJ uses these results to find spatio-temporal deviations of the 
event size distribution from its mean-field form. Results of numerical experiments are 
illustrated in Sect. IVIII In Section IVIIII we further discuss the relation of our results to 
prediction of extreme events. Proofs and necessary technical information are collected in 
Appendices. 
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II. MODEL OUTLINE AND PREDICTION PROBLEM 



Our model combines external driving ultimately responsible for occurrence of events, 
including the extreme ones, a cascade process responsible for redistribution of energy (or 
another appropriate physical quantity such as mass, moment, stress, etc.) within the 
system, and spatial dynamics. We first outline the process of populating a system space Q 
with particles of discrete ranks and then proceed with definition of the observation space 
and events. We assume that Q is an n-dimensional Euclidean space. 

A direct cascade (branching) within a system starts with consecutive injection (immi- 
gration) of particles of the largest possible rank, r max , into the origin 6 fl, which we 
call source. After injection, each particle diffuses freely and independently of the others 
across the space Q. Eventually, it splits into a random number of particles of smaller rank, 
r max — 1, each of which continues to diffuse from the location of the parent and indepen- 
dently of the other particles. These particles split in their turn into even smaller particles, 
and so on. 

At each time instant t > 0, observations can be done on a subspace 7Z t C Q. In this 
paper we assume that lZ t is an affine subspace of dimension d < n. An observed event 
corresponds to an instant when a particle crosses the subspace of observations. Each event 
is characterized by its occurrence time t, spatial location x e TZ t within the observation 
space, and rank r. Model observations at instant t thus consist of a collection of events 
Ct — {ti < t,Xi,Ti), % > 1, referred to as catalog. Extreme event is defined as a sufficiently 
large, although not necessarily the largest, event, r > r , where r is a rank threshold. 

Importantly, the location of TZt within Q is a) not known to the observer, and b) time- 
dependent. One can interpret this as movement of the observation space relative to the 
source, movement of the source relative to the observation space, or combination of the 
two. A principal goal of an observer is to assess the likelihood of the occurrence of an 
extreme event using the catalog C t . It is readily seen that the probability of an extreme 
event increases as the observation space approaches the source and achieves its maximal 
value when the source belongs to the observation space, G TZ t . The distance between the 
observation subspace and the source thus becomes a natural control parameter and allows 
one to reduce the prediction problem to estimating the distance to the source. This latter 
problem is the focus of our study. 

As the observation subspace approaches the source, intensity of the observed events 
increases, larger events become relatively more frequent, clustering and long-range correla- 
tions become more prominent (see Fig. [Hand Sects. lIVllVIIII) . Emergence of these patterns, 
each individually and all together, can be therefore used to forecast an approach of a large 
event; indeed, such a prediction should be understood in a statistical sense. This study is 
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focused on quantitative description of two of these patterns, intensity increase and devia- 
tions from self-similarity, for a classical branching process formally introduced in the next 
section. 

We emphasize that the location and dynamics of the observation space IZt within Q 
depend on details of a particular system of interest and may be hard to estimate or model. 
An important result of this paper is that (i) the information about this unknown dynamics 
can be summarized by a scalar value of the control parameter (distance between the obser- 
vational subspace and the origin); and (ii) knowledge of the control parameter is sufficient 
to solve the prediction problem. 

Finally, it is important to mention that we do not use direct cascade as a dynamical 
model of event formation, which would imply that large events cause smaller ones. We 
merely use this analytically tractable approach to create a hierarchical network of spatially 
distributed particles. A dynamic interpretation of the latter will depend on a concrete 
application, and may include inverse cascading or other physically relevant processes. 

III. MODEL FORMULATION 

We consider an age-dependent multi-type branching diffusion process with immigration 
in R™. The system consists of particles, each of which belongs to a generation k = 0, 1, . . . . 
Particles of zero generation (the largest ones) appear in a system as a result of external 
driving (forcing); we will refer to them as immigrants. Particles of any other generation k > 
are produced as a result of splitting of particles of generation k — 1. Immigrants (k = 0) 
are born at the origin x := (xi, . . . , x n ) = at a constant rate fi; that is, the probability for 
a new immigrant to appear within the time interval of length At is fiAt + o(At) as At — > 0. 
Accordingly, the birth instants form a homogeneous Poisson process with intensity fi. Each 
particle lives for some random time r and then transforms (splits) into a random number 
j3 of particles of the next generation. The probability laws of the lifetime r and branching 
(3 are generation-, time-, and space-independent. We assume that new particles are born 
at the location of their parent at the moment of splitting. 

The particle lifetime has an exponential distribution: 



The conditional probability that a particle transforms into k > new particles (0 means 
that it disappears) given that the transformation took place is denoted by p k . The proba- 
bility generating function (pgf) for the number (3 of new particles is thus 



G(t) := P{r <t} = l-e- xt , A > 0. 



("I 1) 




(III.2) 



k 
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The expected number of offsprings (also called the branching number) is B := E(/3) = h'(l) 



(see e.g. 
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Chapter 1). 

Each particle diffuses in K n independently of other particles. This means that the density 
p(x, y, t) of a particle that was born at instant at point y solves the equation 

!^(e5)^^ iP en,) 

with the initial condition p(x, y, 0) = <5(x — y). The solution of (1III.3I) is given by (81] 

p(x,y,t) = (4 7 r J Dt)-"/ 2 exp|-^^|, |x| 2 = J> 2 . (m.4) 

Accordingly, the density of each particle, given that it is alive at the instant t, is <f>(x.,t) : = 
p(x, 0,t). Naturally, the positions of the particles produced by the same immigrant are 
correlated. This can be reflected by the joint distribution of pairs, triplets, etc. 

The model is specified by the following parameters: immigration intensity fi > 0, branch- 
ing intensity A > 0, diffusion constant D > 0, and branching distribution {pk}, which will 
be often represented by its pgf h(z) or simply by the branching number B. An appropriate 
choice of the temporal and spatial scales allows one to assume \l — 1 and D — 1. 

It is convenient to introduce particle rank r := r max — k for an arbitrary integer r max and 
thus consider particles of ranks r < r max . Particle rank can be considered a logarithmic 
measure of the size. Similar to the analysis of the real-world systems, we sometime only 
consider particles of the first several generations < k < r max — 1, which corresponds to 
the largest ranks 1 < r < r max . Figure CD illustrates the model population. 



IV. SUMMARY OF RESULTS RELATED TO PREDICTION 

We summarize here the study's findings that are most relevant to the prediction problem. 
Recall that the prediction problem consists of assessing the likelihood of an extreme event; 
the latter corresponds to an instant when a sufficiently large particle crosses the observation 
space. The likelihood of an extreme event is thus directly related to the distance between 
the space of observations and the origin. Accordingly, the prediction problem is reduced to 
the estimation of this distance from available data. For that, one should look for increase in 
the intensity of medium-to-large-sized events, as well as upward deviations in the event size 
distribution. We believe that this general idea can be useful in a wide range of models and 
observed systems, not necessarily based on a branching diffusion mechanism. Statistical 
assessment of particular prediction schemes based on this idea is left for a future study. 

All statements below refer to a steady-state of the model (dynamics after a transient). 
All asymptotic statements have been confirmed numerically in finite models. 
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Meanfield self- similarity. Particle ranks, averaged over time and space, have an ex- 
ponential distribution; this is equivalent to a power-law distribution of particle sizes; 
see (jVEED and Fig. El 

Small-size self- similarity. The particle rank distribution at any spatial point is asymp- 
totically exponential as rank decreases, with the exponent index —B; see (1VI.6P and 
Figs. |2] and HI This is equivalent to a power-law distribution of particle sizes with 
power-law index —B. Furthermore, this implies that deviations from self-similarity, 
if any, can be only seen at large ranks (large particle sizes). 

Upward deviations close to the origin. At any point sufficiently close to the origin, 
the particle size distribution deviates from a self-similar power-law form as to have 
a larger number of medium-to-large-sized events. The magnitude of this deviation 
increases with the event size, as well as with dimension of the model space; see (1VI.4I) 
and the upper lines in Figs. [2] and HI 

Downward deviations away from the origin. At any point sufficiently far from the 
origin, the particle size distribution deviates from a self-similar power-law form as 
to have a smaller number of medium-to-large-sized events. The magnitude of this 
deviation increases with the event size and is independent of the model's dimension; 
see (1VI.5|) and the lower lines in Figs. |2] and HI 



5. Exponential decay of event intensity. The intensity of events of any fixed size is 
exponentially decaying away from the origin; see flV.24j) . 

6. Divergence of event intensity at the origin. For models with spatial dimension larger 
than 1, the intensity of sufficiently large events diverges at the origin in a power- law 
fashion; see (TvT^ . dVTM and Fig. E2(b,c,d). 



V. MODEL SOLUTION: MOMENT GENERATING FUNCTIONS 

The model introduced in Sect. IHIl is a superposition of independent branching processes 
generated by individual immigrants. Sections IV Al and IV Bl analyze, respectively, the one- 
point and two-point moments of a particle distribution produced by a single immigrant. 
Then we expand these results to the case of multiple immigrants in Sect. IV CI 
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A. Single immigrant: One-point properties 

1. Moment generating functions 

Let Pk,i{G,y,t) be the conditional probability that at time t > there exist i > 
particles of generation k > within spatial region GcM" given that at time a single 
immigrant was injected at point y. The corresponding moment generating function is 

M k (G, y, t- s) = Y^PkAG, y, t)e s \ (V.l) 

i 

Proposition V.l The moment generating functions M k (G,y,t; s) solve the following re- 
cursive system of non-linear partial differential equations: 

^M k (G, y, t; s) = DA y M k - XM k + A h(M k ^), k > 1, (V.2) 

with initial conditions M k (G, y,0; s) = 1, k > 1, and 

M {G, y, t; s) = (l-P) + Pe s , P := e' xt [ p(x, y, t)dx. (V.3) 

iiTere h(s) is defined by 0111.20 and A y = J2i9 2 /dyf. 
Proof is given in Appendix [A] 



2. The first moment densities 

Let Ak(G,y,t) be the expected number of generation-A; particles at instant t within the 
region G, produced by a single immigrant injected at point y at time t = 0. It is given by 
the following partial derivative (see e.g., 



80|, Chapter 1): 



A k (G,y,t) 



dM k (G,y,t ]S ) 



ds 



s=0 



(V.4) 



Consider also the expectation density A k (x,y,t) that satisfies, for any G C M. n , 



A k (G,y,t) = / A fc (x,y,t)dx. 

' G 



(V.5) 



Corollary V.2 The first moment densities A k (x.,y,t) solve the following recursive system 
of partial differential equations: 



dA k (x,y,t) 
dt 



DA x A k - XA k + XBA k ^ x , k>l, 



(V.6) 
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with the initial conditions Ak(x., y, 0) = 0, k > 1, 

A (x,y,0) = 5(y-x), A (x, y, t) = e~ A V(x, y, t), t > 0. (V.7) 
27ie solution to this system is given by 

A k (x,y,t) = ^i-^ ( x ,y,*) 

- (A ^-" /2 exp(-At-^#l. (V.8) 



k\(AnD) n / 2 * { Wt 



Proof is given in Appendix O It follows from a general result for the higher moments 
obtained in Appendix [Bj 

The system flV.6j) has a transparent intuitive meaning. The rate of change of the ex- 
pectation density A k (x, y, t) is affected by the three processes: diffusion of the existing 
particles of generation k (the first term in the rhs of (IV. 6p ). splitting of the existing parti- 
cles of generation k at the rate A (the second term), and splitting of the generation k — 1 
particles that produce on average B new particles of generation k (the third term). 

To obtain the solution for the entire population, we sum up the contributions from all 
generations: 

°° „-AI(l-B) / |„|2 \ 

A(x, 0, t) = £ A fc (x, 0, t) = e-^~ B ) p(x, 0, t) = exp -±±- . (V.9) 

This formula emphasizes the role of the branching parameter B: in subcritical case, B < 1, 
the population extincts exponentially; in supercritical case, 5 > 1, the population grows 
exponentially; in critical case, B — 1, the expected number of particles remains the same 
(steady state) and is given by the diffusion density p(x, 0, t). 



B. Single immigrant: Two-point properties 

1. Moment generating functions 

Let pk l! k 2 ,i,j{Gi,G2,y,t) be the conditional probability that at instant t > there exist 
i > particles of generation k\ > within region Gi C M n and j > particles of generation 
A>2 > within region G2 C M n given that at time a single immigrant was injected at point 
y. Assume that G\ and G2 do not overlap. The corresponding moment generating function 
is 

M kuk2 (G 1 ,G 2 ,y,t;s 1 ,s 2 ) = ^ V^m^A^^ G 2 ,y,t)e isi+jS2 . (V.10) 

i,j>0 
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Proposition V. 3 The moment generating functions M kltk2 (G\, G 2 , y, t; s±, s 2 ) solve the 
following recursive system of non-linear partial differential equations: 

d 

— M klM = DA y M klM - \M klM + A h(M kl ^ 1M ^), k u k 2 > 1, (V.ll) 

with the initial conditions 

M klM (G 1 ,G 2 ,y,0;s 1 ,s 2 ) = l, h,k 2 >l, (V.12) 

M ,o (G u G 2 , y, t; s u s 2 ) = P^ + P 2 e s * +1-P 1 -P 2 , (V.13) 
M 0tk (G 1 , G 2 , y, t; Sl , s 2 ) = {M k (G 2 , y, t; s 2 ) - e~ xt ) + (e- xt - P 1 ) + P ie s \ (V.14) 

where Pi := e~ xt f G jo(x, y, t)dx, i = 1,2. Here, as before, h(s) is defined by 2 j) and 
Proof is given in Appendix [D] 



2. Moments 



Consider the expected value A kltk2 (Gi,G 2 ,y,t) of the product of the number of 
generation-/^ particles in region G\ and number of generation-/c 2 particles in region G 2 
at instant t, produced by a single immigrant injected at point y at time t = 0. It is given 
by the following partial derivative 



A kuk2 (G 1 ,G 2 ,y,t) :-- 



d 2 M klM (Gi, G 2 , y, t; si, s 2 ) 



ds\ds 2 



(V.15) 



Sl=S2=0 



We notice that the expectations A kl (Gi,y,t) and A k2 (G 2 ,y, t) of (1V.4|) can be repre- 
sented as 



A kl {G u y,t) :-- 



dM klM (Gi, G 2 , y, t; si, s 2 ) 



and 



A k2 (G 2 ,y,t) :-- 



dsi 



dM klM (Gi, G 2 , y, t; si, s 2 ) 



(V.16) 



Sl=S2=0 



OS-} 



(V.17) 



=s 2 =0 



Consider also the expectation density A^^Xi, x 2 , y, t) that satisfies, for any nonoverlap- 
ping G 1 ,G 2 cW n 



A klik2 (G u G 2 ,y,t) = A kl > fc 2 (x 1 , x 2 , y, t)<ix 1 (ix2. 

>G 2 JG! 



(V.18) 
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Corollary V.4 The moment densities A^^ = ^4fci,fc 2 ( x i> x 2, y, t) solve the following re- 
cursive system of partial differential equations: 

— 2 = -D AyAfafa — A Ak 1: k 2 + A B Afa-ifa-i 

+ A /*"(!) ^(x^^^Xa), h,k 2 >l, (V.19) 



iint/i £/ie initial conditions 



A fcl)fe ( Xl ,x 2 ,y,0) = 0, fei, fc 2 > 1, (V.20) 
A),fc(xi, x 2 , y, i) = 0, k > 0, t > 0, (V.21) 



and v4fc(x) = Afe(x, y, t) given by (IV. 8j) . 
Proof is given in Appendix [Ej 

C. Multiple immigrants 

Here we expand the results of the Sect. IV Al to the case of multiple immigrants that 
appear at the origin according to a homogeneous Poisson process with intensity fi. The 
expectation Ak of the number of particles of generation k is given by 

A(x,t) = / A k (-x, 0, s)/ids 



o 



ii (A B) k r* sk _ n/2 exp r _ A s _ jx^_ I ^ 

Jo I 4L>sJ 



fc! (4 7rD) n/2 

The steady-state spatial distribution corresponds to the limit t — > oo: 

A(x) := A.(x,oo) = , 2 ^ (Ai?) l ( J^rV /2 /C f |x| I . (Y.23) 



A;!(4 7r J D) ri/2 V4^A; ^ 1 V £ 

Here z/ = & — n/2 + 1 and lf„ is the modified Bessel function of the second kind (see 
Appendix |G|) . Introducing the normalized distance from the origin z := |x| ^f\jD we 
obtain 

fj, (B\ k {2nD^ ~ n/2 



*M=mi — 1 2 " /cw ' (v - 24) 



For odd n, there are explicit expressions for K„(z) (Appendix |G| Eqs. ( 1G.2j) . (1G.3j) ). In 
particular, we have 



Ao(z) = -j==e-\ for n = 1, A (z) = J A -^-e"*, for n = 3. (V.25) 
V4 D A V 4 7r z 
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From (1V.24|) and the asymptotic behavior of K v (z) as z — > (Appendix [GJ Eq. (1G.5j) ) it 



follows that 

J oo, for v < 0, i.e., k < n/2 - 1 
hm A(^) = < , (V.26) 

z^o I const < oo, for z/ > (J, i.e., fc > n/2 — 1. 

Thus, in a model with spatial dimension n > 2, the elements of several lowest generations 
(k < n/2 — 1) have an infinite concentration at the origin. 



D. Alternative model representation 

In this section we derive a system of equations for the steady-state expectations Ak{x) 
using the radial symmetry of the problem. By integrating the equation flV.6j) from t = 
to oo, we obtain 

DA x A(x) - AA(x) + ABA-i(x) = 
since A&(x, y, oo) = 0. We now rewrite this equation in terms of the normalized distance 



from the origin, z := \x\y/X/D, using the fact that -4fc(x) = Ak(z) as soon as |x| 



\z\ 



Kit) + ^ A' k (z) - A k {z) + BAk-i(z) = 0. (V.27) 



We notice, furthermore, that one can rewrite the expectation densities ( IV. 81) as a function 
of z, which results in Ak(z) = Afe(x, 0, t). It is then readily seen that 

Ak(z) = ~2T z Ak ~ l{z) - (v ' 28) 

The same recursive system holds for Ak(z), which is shown by integrating the last equation 
with respect to time. 



VI. PARTICLE RANK DISTRIBUTION 

We analyze here the particle rank distribution; recall that the rank is defined as r = 
r max — k, where k is the particle's generation. A self-similar branching mechanism that 
governs our model suggests an exponential distribution of particle ranks. Indeed, the 
spatially averaged steady-state rank distribution is a pure exponential law with index B: 



p poo 

Ak : = / / Afc(x, 0, t)ixdt dx 

JR n JO 

t 10 I f\ 4\h -At jj. P Tjk 



n (Xt) k e- Xt dt= ^-B k <xB- r . (VI.l) 
k\ J a 
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Remark VI. 1 Our use of the term "self- similar" with respect to the exponential distri- 
bution, often seen in physical literature, requires some explanations. As we mentioned 
earlier, the particle rank serves as a logarithmic measure of its size. Thus, the exponen- 
tial distribution of ranks corresponds to the power law distribution of sizes; hence the term 
"self-similarity" . 

To analyze rank- and space-dependent deviations from the pure exponential distribution, we 
will consider the ratio 7&(x) between the number of particles of two consecutive generations: 

For the purely exponential rank distribution, Ak(x) = cB k , the value of 7fc(x) = 1/B is 
independent of k and x; while deviations from the pure exponential distribution will cause 
7^ to vary as a function of k and/or x. Plugging flV.24j) into flVI.2j) we find 

2(fc + l) K v (z) 
7t(X>= B z K^SJH ' (VL3) 



where, as before, z := |x| y/\/D and v = k — n/2 + 1. 

Proposition VI. 2 The asymptotic behavior of the function jk( z ) is given by 

v < 0, 

lim 7fc (z) = { 1 M n ( VL4 ) 




7^(2) ~ — — , z — > oo, fixed k, (VI. 5) 

1 / n \ 

lk{z) ~ — ( 1 H ), k oo, fixed z. (VI. 6) 



Proof and explicit rates of divergence in (1VI.4|) are given in Appendix [Fl 

Proposition IVI.2I describes the spatio-temporal deviations of the particle rank distri- 
bution from the pure exponential law (IVI.lj) . We interpret below each of the equations 
(1VI.4I) - (1VI.6|) in some detail. Eq. f ]VI.6j) implies that at any spatial point, the distribu- 
tion asymptotically approaches the exponential form as generation k increases (rank r 
decreases). In other words, the distribution of small ranks (large generation numbers) is 
close to the exponential with index —B; thus the deviations can only be observed at the 
largest ranks (small generation numbers). Analysis of the large-rank distribution is done 
using Eqs. ( 1VI.4I) and (I VI. 51) . Near the origin, where the immigrants enter the system, 
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Eq. (1VI.4|) implies that 7&(z) > jk+i( z ) > 1/-B for z/ > 0. Hence, one observes the upward 
deviations from the pure exponential distribution: for the same number of rank r particles, 
the number of rank r + 1 particles is larger than predicted by the exponential law. The 
same behavior is in fact observed for v < (see Appendix [Fj Eq. flF.5j) ). In addition, 
for v < the ratios 7&(V) do not merely deviate from 1/B, but diverge to infinity at the 
origin. Away from the origin, according to Eq. (IVI.5j) . we have jk( z ) < Jk+i( z ) < 
which implies downward deviations from the pure exponent: for the same number of rank 
r particles, the number of rank r + 1 particles is smaller than predicted by the exponential 
law. 

Figure [2] illustrates the above findings; it shows the distribution of particles for the 
largest ranks at different distances from the origin. One can clearly see the transition from 
downward to upward deviation of the rank distributions from the pure exponential form 
as we approach the origin. Notably, the magnitude of the upward deviation close to the 
origin (the upper line in all panels) strongly increases with the model dimension n. 

VII. NUMERICAL ANALYSIS 

Our analytical results and asymptotics are closely reproduced in numerical experiments 
with finite number of generations, limited spatial extent, and spatial averaging (unavoidable 
when working with observations). Here, to mimic the ensemble averaging, the numerical re- 
sults have been averaged over 4000 independent realizations of a 3D model with parameters 
fj, = \ = 1, D = 1, and B = 2. 

First, we check the exponential rank distribution of (IVI.lj) . Figure [3] shows the observed 
spatially averaged particle rank distribution. The exponential form (1VI.1I) is indeed well 
reproduced. 

Next, we see how the spatial averaging affects the rank distribution. Figure H] shows 
the rank distribution at t = 30 at various distances to the origin. The spatial averaging 
has been done within spherical shells (space between two concentric spheres) of a constant 
volume V = 5. Thus, here we see an observable counterpart of the theoretical distributions 
shown in Fig. [2b. Although the spatial averaging somewhat tapers off the upward bend at 
the largest ranks close to the origin, the predicted transition from the downward to upward 
bend is clearly seen. 

Figure |5] illustrates in more detail how the spatial averaging affects the upward bend 
in a 3D model. It shows the particle rank distributions at t = 30 spatially averaged over 
spheres of different volumes centered at the origin. The upward bend is prominent for the 
spheres with volumes V < 5; and it gradually disappears within larger spheres in favor of 
an exponential distribution observed after a complete spatial averaging. Notably, the pure 
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exponential distribution can be only achieved by averaging over all events in the model 
(V= oo). 



VIII. DISCUSSION 

This work is motivated by the problem of prediction of extreme events in complex 



systems. Our point of departure is the four types of premonitory patterns [15J], previously 
found in models and observations. We propose here a simple mechanism and a single 
control parameter for all these patterns. 

Quantitative analysis is performed here for a classical model of spatially distributed pop- 
ulation of particles of different sizes governed by direct cascade of branching and external 
driving (see Sect. IIIII) . In the probability theory this model is known as the age-dependent 
multitype branching diffusion process with immigration [sQ|. We consider here a new scope 
of problems for this model. We assume that observations (detection of particles) are only 
possible on a subspace of the system space while the source of external driving (origin) re- 
mains unobservable, as is the case in many real-world systems. The natural question under 
this approach is the dependence of the process statistics on the distance to the source. A 
complete analytical solution to this problem, in terms of the moments with respect to the 
particle density, is given by Proposition IV. II In addition, the correlation structure of the 
particle field can be found using Proposition IV.3I 

It is natural to consider rank as a logarithmic measure of the particle size. The ex- 
ponential rank distribution derived in Eq. fIVI.lj) corresponds to a self-similar, power-law 
distribution of particle sizes, characteristic for many complex systems. The self-similarity 
in our model, as well as in the real-world systems, is only observed after global spatial 
averaging in a steady-state. Proposition IVI.2I and Fig. [2] describe space-dependent devia- 
tions from the self-similarity. Recall that an extreme event in our system is defined as an 
observation of a particle of sufficiently large size. As the source approaches the observa- 
tion subspace, the probability of an extreme event increases. Our results are thus directly 
connected to prediction: When the location of the source changes in time approaching the 
subspace of observation (or vice versa), the increase of event intensity and the downward 
bend in the event size distribution becomes premonitory to an extreme event. The numer- 
ical experiments confirm the validity of our analytical results and asymptotics in a finite 
model. 

Our model exhibits very rich and intriguing premonitory behavior. Figure [1] shows 
several 2D snapshots of a 3D model at different distances from the source. One can see that, 
as the source approaches, the following changes in the background activity emerge: a) The 
intensity (total number of particles) increases; b) Particles of larger size become relatively 
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more numerous; c) Particle clustering becomes more prominent; d) The correlation radius 
increases. All these premonitory changes have been independently observed in natural 
and socioeconomic systems. Here they are all determined by a single control parameter - 
distance between the source and the observation space. 



The abovementioned premonitory patterns closely resemble universa' 
els of statistical physics in a vicinity of second order phase transition 



properties of mod- 
84|, percolation 



82 



models near the percolation threshold [85l. 186], or random graphs prior to the emergence of 



a giant cluster |87H89j. In these models, the approach of an extreme event, usually referred 



to as critical point, and the emergence of premonitory patterns, called critical phenomena, 
correspond to an instant when a control parameter crosses its critical value. In statistical 
physics a typical control parameter is temperature or magnetization; in percolation it is 
the site or bond occupation density; in a random graph — the probability for two vertices 
to be connected. The theory of critical phenomena [83j quantifies system's behavior at the 
critical value of the corresponding control parameter. The remarkable power of this theory 
is connected to the fact that very different systems demonstrate similar behavior near to 
criticality. More precisely, when the control parameter is close to its critical value, the 
system sticks to one of just a few types of possible limit behaviors, each being described 
by an appropriate scale-invariant statistical field theory. In particular, each limit behavior 
corresponds to the asymptotic power-law size distribution of system observables with a 
characteristic value of critical exponent. 

We focus here on a problem inverse to that considered by the critical phenomena theory: 
Estimating the deviation of a control parameter from the critical value using the observed 
system behavior. The motivation for this is coming from environmental, geophysical, and 
other applied fields where one faces a problem of assessing the likelihood of occurrence of 
an extreme event associated with a critical point. We formulate and solve such a prediction 
problem for a spatially embedded cascade process, which enjoys both the mean-field self- 
similarity and realistic premonitory time- and space-dependent deviations from the latter. 
The methods developed in this paper may provide a framework for studying predictability 
of extreme events in complex systems of arbitrary nature. 
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Appendix A: Proof of Proposition IV. II 



We will need the following calculus lemma that is readily proven by using the definition 
of derivative: 

Lemma A.l Let f(z), g(z), z G R be continuous functions such that the definite integral 
G(t) = j f(z) g(t — z)dz exists. We also assume that g{z) is differentiable. Then, 

jG{t) = J*f(z)g'(t-z)dz + f(t)g(0). 



There are two possible scenarios for the model development up to time t. In the first 
one, the initial immigrant will not split; the probability for this is P = e~ xt . In the second 
one, the initial immigrant will split at instant < u < t; the probability of the first split 
within the time interval [u, u + du] is Xe~ xt du + o(du) as du — > 0. The spatial position of the 
split is given by the diffusion density p(x, y,u). If the immigrant splits, the composition 
property of generating functions gives M k = h[M k _ 1 }. Integrating over all possible split 
instants and locations, we obtain 

M k (G,y,t;s) = e- xt + [ dy'l du Ae" A >(y', y, u) h[M k ^(G, y', t - u; s)}. (A.l) 

JR n JO 

Here the first and the second terms correspond to the first and second scenarios, respec- 
tively. Using the new integration variable z = t — u, we write 

M k (G, y, t; s) = e~ xt + e~ xt [ dy' f du Xe x ^p(y', y, u) h[M k ^(G, y', t - u; s)} 

JR n JO 

(l+ [ dy' [ dz\e Xz p(y',y,t-z)h{M k _ 1 (G,y',z;s)]). 
\ Jr" Jo J 



e~ xt 



Now we take the derivative with respect to t of both sides and apply Lemma A.l using 
the fact that p(y', y, 0) = 5(y' - y) and (d/dt - DA y )p = : 

^M fc (G,y,t;s) = —\M k (G,y, t; s) 
+e- xt \ [ dy' [ dz Xe Xz h[M k ^(G, y', z; s)]DA y p(y', y,t-z) + Xe xt h[M k ^(G, y, t; s)] 

jR n Jo 

Taking the operator A y out of the integration signs, we find 

^- f M k (G, y, t; s) = DA y M k - XM k + X h[M k ^\. 
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It is left to establish the initial conditions. Since we start the model with a particle of 
generation k = and the distribution of splitting is continuous, at t — there are no other 
particles with probability 1. Hence, Mk(G,y, 0; s) = 1 for all k > 1. For generation k = 0, 
we can only have one or no particles at time t > 0. The probability to have one particle 
is given by the product of probabilities that there was no split up to time t and that the 
particle happens to be within region G at time t: P = e~ xt j G p(x, 0, £)<ix. The probability 
to have no particles is then (1 — P). This implies f]V.3|) . □ 



Appendix B: Moments in one-point system 



For any natural number j, consider the j-th moment A\/\G, y,t) of the number of 
generation-fc particles at instant t within the region G, produced by a single immigrant 
injected at point y at time t = 0. It is given by the following partial derivative (see e.g., 



80|, Chapter 1): 



^\G,y,t):- 



&M k {G,y,t; s) 



dsi 



s=0 



(b.i; 



Corollary B.I The moments A^\G,y,t) solve the following recursive system of partial 
differential equations: 



^\G,y,t) = DA y A^-\A^ + \ 



E 



mi\m2>. . . .nij 



-/ m) (l)fl 



Ad) 
A k-1 



i=l 



where m = m\ + ■ ■ ■ + rrij and the sum is over all partitions of j , i.e., values of mi 
such that mi + 2m2 + • • • + jmj = j , with the initial conditions 

Af(G, y,0) = 0,*>1, 



4 j) (G,y,0) = / 5(y-x)rfx, 
Jg 

A^{G,y,t) = e~ xt [ p(x,y,t)dx, t > 0, 



(B.2) 



,...,mj 



(B.3) 
(B.4) 

(B.5) 



and 



G 



s=l 



E 



711 



n — i) 



Pn- 



Proof: The validity of (IB. 2 1) follows from Proposition IV. II Namely, applying the operator 
&* /dsi (-)\s=o to both sides of ( IV. 21) . changing the order of differentiation, and using Faa di 
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Bruno's formula for the j-th derivative of a composition function, one finds, for each k > 1, 



dM k (G,y,t;s) 
dt 



d' 



s=0 



dsi 



- [DA y M k - XM k + A h(M k ^)} 



s=0, 



d 


~&M k (G,y,t;s) 






dt 




s=0- 





r , x d j M k d j M k x # , 
9s-? <9s-? <9s^ 



s=0 



|4 B (G,y,«) 



+ ME 



DAy—— - X— - 

y dsi dsi 



m 1 !m 2 ! . . .rrijl 



3 ( M (i) 

i=l 



DA y Af - XA^> + A 



rO) 



E 



s=0 

(j) \ m i 



m 1 !m 2 ! . . .rrij 



J I AW 
-h^ (1)11 1 



, mj 



where m = m! + • ■ ■ + rrij and the sum is over all partitions of j, i.e., values of mi, 
such that m\ + 2m 2 + ■ ■ • + jwij = j. The initial conditions are established by applying 
the operator $ f ds* (-)\ a= Q to both sides of (1V.3jl and using the definition of A^'(G,y,t) in 

(EU). 



Appendix C: Proof of Corollary IV. 21 



For j = 1 , the equation in Corollary IB. II simplifies to 





dt 



A k (G,y,t) = DA y A k - XA k + XBA k _ x . 



Using the definition of Afc(x, y, t) given in (IV. 5j) . one obtains for each k > 1 



d 

— A fe (x, y, t) = DAyAfc - AA fc + ASA^i. 

It is left to use the translation property A k (x., y, t) = A k (x. — y, 0, t) to change A y to A x . 
The validity of general solution (lV.8j) is proven by induction using the fact that 



d_ 

dt 



DA Y + A 



A = 0. 



The last equality in (PvT8]) follows from flB~5]) and fllTLil) . 



□ 
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Appendix D: Proof of Proposition IV. 31 



The proof of Proposition IV. 31 follows the line of the proof of Proposition IV.ll There 
are two possible scenarios for the model development up to time t. In the first one, the 
initial immigrant will not split; the probability for this is P = e~ xt . In the second one, the 
initial immigrant will split at instant < u < t; the probability of the first split within 
the time interval [u, u + du] is \e~ xt du + o(du) as du — > 0. The spatial position of the 
split is given by the diffusion density p(x, y, u). If the immigrant splits, the composition 
property of generating functions gives M kltk2 — M^fci-i,fe 2 -i]- Integrating over all possible 
split instants and locations, we obtain 

M klM (Gi,G 2 ,y,t;s 1 ,s 2 ) 

= e~ xt + [ dy' [ du Ae- A >(y / , y, u) h [M fel _ ljfe2 _i G 2, y', t - u; s u s 2 )\ . 

JR n JO 

Here the first and the second terms correspond to the first and second scenarios, respec- 
tively. Using the new integration variable z = t — u, we write 

M klM (G 1 ,G 2 ,y,t;s 1 ,s 2 ) 

= e~ xt + e- xt [ dy' [ du Xe x ^p(y', y, it) h [M kl ^ k2 ^ (G u G 2 , y', t - it; Sl , s 2 )] 

JR n JO 



e -xt 



(l+ [ dy' [ dz\e Xz p(y',y,t- z)h[M kl _ 1M _ l (G 1 ,G 2 ,y',z;s u s 2 )}) 
\ Jr" Jo J 



Now we take the derivative with respect to t of both sides and apply Lemma [A. II using the 
fact that p{y', y, 0) = S(y' - y) and (d/dt - DA y )p = 0: 



d 

—M klM (G 1 ,G 2 ,y,t;s 1 ,s 2 ) = ~\M klM {G x , G 2 , y, t; s u s 2 ) 



+ e~ xt 



[ dy' [ dz Xe Xz h [Mfc-i,*,,-! (G X ,G 2 , y', z; s u s 2 )} DA y p (y', y,t - z) 
yr Jo 

+ Xe xt h [Affe-i^-i (G h G 2 , y, *; s x , s 2 )\ 

Taking the operator A y out of the integration signs, we find 
d 

—M klM (Gi, G 2 , y, t] si, s 2 ) = DA y M klM - \M klM + X h[M kl -i,k 2 -i\- 

It is left to establish the initial conditions. Since we start the model with a particle of 
generation k = and the distribution of splitting is continuous, at t — there are no other 
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particles with probability 1. Hence, M kuk2 {pit G 2 , y> 0; s ii s 2) = 1 for all ki,k 2 > 1. For 
generation k\ = k 2 = 0, we have three possibilities: the initial immigrant has not split and 
is in G\ (i — l,j — 0), the initial immigrant has not split and is in G 2 (i = 0,j = 1), 
and neither (i = 0,j = 0), with corresponding probabilities of Pi, P 2 , and 1 — P\ — P 2 , 
respectively. This implies (1V.13|) . 

For generation k\ = and k 2 — k > 1, we again have three possibilities: the initial 
immigrant has not split and is in G\ (i — l,j — 0), the initial immigrant has not split 
and is not in G\ (i = 0, j = 0), and the initial immigrant has split (i = 0,j > 0), with 
corresponding probabilities of Pi, e~ xt — Pi, and 1 — e~ xt , respectively In the last case, the 
number of the 0-th generation particles in Gi is with probability 1 while the information 
on the k-th generation particles in G 2 is given by 

/ dy' f du\e- Xu P (y',y,u)h[M k _i(G 2 ,y',t-u;s 2 )}. 

JR n JO 

From flA.ll) . we see that the above expression equals Mk{G 2 ,y,t; s 2 ) — e~ xt . This implies 
fTVUil) . We notice that setting s 2 = in fTVU3|) and fTVUil) each yields (1 - Pi) + Pie Sl as 
it should (cf. (TP]) ). □ 



Appendix E: Proof of Corollary IV. 41 



The validity of (IV. 19[) follows from Proposition IV. 31 and the definition of 
4fci,*a(Gi> G 2, y, t), A fel)fe2 (x 1; x 2 , y, t). Formally, applying the operator d 2 /dsids 2 (-)\ Sl=S2=0 
to both sides of ( IV. 1 If) and changing the order of differentiation, one finds, for each 
h,k 2 > 1, 



2 



ds 1 ds 2 



dM, 



ki,k 2 



dt 



d 2 



Sl = 



d 


~d 2 M klM 






dt 


dsids 2 


si=S2=0 





DA 



dsids 2 
d 2 M kl 



[DA y M klM - XM klM + A/i(M fel _ 1)fe2 _!)] 



sj=s 2 =0 



k 2 



A 



d 2 M klM 



+\h"(Mi 



y dsids 2 dsids 2 

i,fc 2 -i 



+ Xfi (M kl _i M _i 



d 2 M kl . 



dsids 2 



k\ — l,k2 — l, 



dsi 



ds 5 



31=32=0 



9 -r 

— = DA y A kuk2 



\A klM + APA fcl _i ifc2 _i + Xh" (l)A kl ^i(Gi)A k2 ^i(G 2 ). 



The system (1V.19|) readily follows now from the definition of ^^(xi, x 2 , y, t). The initial 
conditions (lV.20p - flV.21|) are established by applying the operator d 2 / dsids 2 (-)\ sl=S2=0 to 
both sides of (IV. 121) - (IV. 14^ and using again the definition of ^^^(xi, x 2 , y, t). □ 
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Appendix F: Proof of Proposition IVI.2I 



The asymptotic (IVI.5I) readily follows from flG.4j) . To prove flVI.6j) . let r u (z) 
K v {z)/ K v+ x(z). From (IG.ip one finds that 



and furthermore 



K u+1 {z) 
K(z) 

z 1 



K v - X {z) 2u 
KJz) z 



(F.l) 



(F.2) 



2 ur v (z) 2v 

From monotonicity of K u (z) with respect to the index v > it follows that r u (z) < 1 for 
v > 0. Accordingly, the first term in the rhs of (IF.2I) goes to zero as k — > oo. Hence, 

1 



z 



lim 

fc->oo 2 v r,,(z) 



1. 



or r\, m ~ 



2z/ 



— )■ oo. 



(F.3) 



To complete the proof of (1VI.6j) . we use this asymptotic in f ]VI.3j) . Finally, we prove (1VI.4j) . 
In fact, we will derive a stronger result showing the asymptotics of r„(z) and 7„(z) as z — > 0. 
To find the asymptotics for r l/ (z), we use (IG.5P for all possible combinations of signs for v 
and v + 1. We take into account that by definition v can only take values {i, i + l/2}, e z. 



r v (z) 



K v {z) 
K v+1 {z) 
K-i(z) 



T(-u) ( 2 \-v-(-v-i) 



Mi) 



K (z) 

1/2 



K. 



«l/ 2 (*) 

^o(^) 

K!{z) 

K v {z) 



\ K v+1 {z) 

Combining this with (IVI.3j) we find 



r(-iz-i) 
[z(H2/z)-j)]^ 



z(ln(2/z)- 7 ) 
ro+i) \z) 



2(-v-l)/z, v < -3/2, 
— {z In z)^ 1 , f = —1, 
1, i/ = -1/2, 



(F.4) 



— z In z, 
z/{2u), 



f 4 

Bz 2 



lu{z) 



2(k+ 1) 
Bz 



rJz) ~ < 



(z/ + n/2) (-1/ 

("-2) 
Bz 2 In z' 
n-1 
82 ' 
n In g 

- fl + — ) 



v = 0, 
1/ > 0. 



f < -3/2, 

^ = -1/2, 
z/ = 0, 
1/ > 0. 



(F.5) 



One can see that for v < the ratio j v (z) diverges at the origin. The rate of divergence 
increases monotonously from In z to z~ 2 with the absolute value of v. 

Appendix G: Properties of K v 



Here we summarize some essential facts about the modified Bessel function of the second 
kind K v (z). The sources of this as well as further information about K v (z) are handbooks 
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90|, Chapters 9, 10 and [9l|, Sect. 8.4. The function K v can be defined as a decreasing 



solution of the modified Bessel differential equation 



x 



2 y" + xy' - (x 2 + v 2 ) y = 0. 



The function K v {z) exponentially decreases as z — > oo and diverges at z — 0. In 
addition, K- V {z) = K v (z) and 



K u+1 (z) = K u ^(z) + —K u (z). (G.l) 

z 



For integer k > we have 



V 22; ml [k — my. (2z) m 



and in particular 



K±i/2(z) = Jf z e-; # 3/2 (*) = e-. (G.3) 



For arbitrary fixed v and 2 3> v 



K u (z)~J^- e - z , z^oo. (G.4) 



The asymptotic behavior at z = is given by 







(!) ■ 


\v\ 







where 7 f» 0.577 is the Euler-Mascheroni constant. 



log ( - ) - 7, 1/ = 0, 
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FIG. 1: Example of a 3D model population. Different panels show 2D subspaces of the model 3D 
space at different distances |x| to the origin. Model parameters are fi = A = 1, D = 1, B = 2. 
Circle size is proportional to the particle rank. Different shades correspond to populations from 
different immigrants; the descendants of earlier immigrants have lighter shade. The clustering of 
particles is explained by the splitting histories. Note that, as the origin approaches, the particle 
activity significantly changes, indicating the increased probability of an extreme event. 
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FIG. 2: Expected number Ak(z) of generation k particles at distance z from the origin (cf. 
Proposition IVI.2"]) , The distance z is increasing (from top to bottom line in each panel) as z = 
1CT 3 , 2, 5, 10, 20. Model dimension is n = 1 (panel A), n = 3 (panel B), n = 5 (panel C), and 
n = 10 (panel D). Other model parameters: \i = A = 1, D = 1, B = 2, r max = 21. One can 
clearly see the transition from downward to upward deviation of the rank distributions from the 
pure exponential form as we approach the origin. Notably, the magnitude of the upward deviation 
close to the origin (the upper line in all panels) strongly increases with the model dimension n. 
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FIG. 3: Spatially averaged particle rank distribution at t = 30. The distribution is averaged 
over 4000 independent realizations of a 3D model with parameters fi = X = 1, D = 1, B = 2, 
r max = 10. One can clearly see the exponential rank distribution of Eq. (IVI.lj) . 
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FIG. 4: Particle rank distribution at t = 30 and fixed distance z from the origin (cf. Proposition 
IVI.2|h The distribution is averaged over 4000 independent realizations of a 3D model with pa- 
rameters (j, = A = 1, D = 1, B = 2, r max = 10. Different lines correspond to different distances 
(from top to bottom): z = 0, 2, 4, 6, 8. Spatial averaging is done within spherical shells of constant 
volume V = 5 with inner radius z. One can clearly see that the rank distribution deviates from 
the pure exponential form, which corresponds to a straight line in the semilogarithmic scale used 
here. One observes downward deviations at large distances from the origin, and upward deviations 
close to the origin. 
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FIG. 5: Particle rank distribution at t = 30 in a 3D model. The distribution is spa- 
tially averaged over spheres of volume V centered at the origin with (from top to bottom): 
V = oo,500,100,200,5,l,0.01. Model parameters: fj, = A = 1, D = 1, B = 2, r max = 10. The 
upward deviations from the exponential distribution (a straight line) are fading away with the 
extent of the spatial averaging. 
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